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Abstract 

The original contributions of this paper are twofold: a new understanding of the influence of noise 
on the eigenvectors of the graph Laplacian of a set of image patches, and an algorithm to estimate a 
denoised set of patches from a noisy image. The algorithm relies on the following two observations: 

(1) the low-index eigenvectors of the diffusion, or graph Laplacian, operators are very robust to 
random perturbations of the weights and random changes in the connections of the patch-graph; and 

(2) patches extracted from smooth regions of the image are organized along smooth low-dimensional 
structures in the patch-set, and therefore can be reconstructed with few eigenvectors. Experiments 
demonstrate that our denoising algorithm outperforms the denoising gold- standards. 
Keywords: graph Laplacian, eigenvector perturbation, image patch, image denoising 



1. Introduction 

1.1. Problem statement and motivation 

The flrst goal of this work is to study experimentally the perturbation of the eigenvectors of the 
graph Laplacian constructed from a set of image patches. The second goal of the paper is to devise 
a novel method to jointly denoise an image and estimate the eigenvectors of the patch- graph of the 
clean image. 

Recent work in computer vision and image processing indicates that the elusive quest for the 
"universal" transform has been replaced by a fresh new perspective. Indeed, researchers have re- 
cently proposed to represent images as "collage" of small patches. The patches can be shaped into 
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square blocks [1] or into optimized contours that mimic the parts of a jigsaw puzzle [2]. These 
patch-based appearance models are generative statistical models that can be learned from images, 
and are then used to represent images. All the patch-based methods [e.g., 3, 4, 5, 6, 7, 8, 9, 10, 
and references therein] take advantage of the following fact: the dataset of patches, whether it is 
aggregated from a single image, or a library of images, is a smooth set. A consequence of this ob- 
servation is that the eigenvectors of the graph Laplacian [11, 12] provide a natural basis to expand 
smooth functions defined on the set of patches. The image intensity - seen as a function on the 
patch-set - is obviously a smooth function. As a result, the eigenvectors of the graph Laplacian 
yield an efficient (measured in terms of sparsity, for instance) representation of the original image 
from which the patch-set is computed. 

In practice, images are often corrupted by noise, and the eigenvectors are to be computed from a 
graph of noisy patches. The problem becomes: what is the effect of the noise on the perturbation of 
the eigenvectors? Theoretical results [13] provide upper bounds on the angle between the original 
and the perturbed eigenspaces. Unfortunately, as noted by Yan et al. [14], in the context of 
the perturbation of the first non trivial eigenvector of the graph Laplacian, these bounds usually 
overestimate the actual perturbation. In addition, the bounds depend on the separation of the 
eigenvalues, a quantity that is difficult to predict for image patches. Finally, bounds on the angle 
between invariant subspaces cannot be readily translated in terms of the effect of the perturbations 
on the geometric features encoded by the eigenvectors, or the ability of the perturbed eigenvectors 
to approximate the original image. 

1.2. Outline of our approach and results 

The original contributions of this paper are twofold: a new understanding of the influence of 
noise on the eigenvectors of the graph Laplacian and a new method to denoise image patches. 

Our approach relies on a series of experiments aimed at understanding the effect of noise on 
the perturbations of the eigenvectors of the patch-graph. In particular, we show that the low 
order eigenvectors of the graph Laplacian are robust to random variations in the geometry of the 
graph and in the edge weights. These results are connected to similar results in spectral geometry 
quantifying the perturbations of the eigenfunctions of the Laplace-Beltrami operator as a function 
of changes in the metric of the manifold [15, 16]. Equipped with this experimental understanding 
of the robustness of the low index eigenvectors of the graph Laplacian, we propose an iterative 
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procedure to jointly estimate the eigenvectors and the original image. We evaluate our novel 
approach with numerous experiments on a variety of images, and we compare our approach to 
some of the denoising gold-standards. Our approach systematically outperforms these algorithms. 



2. The manifold of patches from a single image 



We consider an image u{x), of size N x N. We extend the image by even symmetry when the 
pixel location x = (i, j) becomes close to the border of the image. We first define the notion of a 
patch. 

Definition 1. Let x^ = {hj) be a pixel with linear index n = i x N + j . We extract an m x m 
block, centered about Xn , 



u{i — j — m/2) ••• u{i — j + m/2) 



(1) 



u{i + m/2^j — m/2) ••• u{i + m/2^j + m/2) 

where, without loss of generality, we take m to he an odd integer, and m/2 is the result of the 
Euclidean division of m by 2. By concatenating the columns, we identify the mxm matrix (1) with 
a vector in W^'^ , and we define the patch u{xn) by 



u{Xn) 



Ui{Xn) 









u{i — m/2, j — m/2) 



u(i + m/2, j + m/2) 



(2) 



As we collect all the patches, we form the patch-set in R"^ . 

Definition 2. The patch-set is defined as the set of patches extracted from the image u, 

V = {u{xn),n = l,2,...,N^}. 



(3) 



Before we start exploring the effect of noise on the patch-set 'P, let us pause for an instant and 
ask ourselves if P is a smooth set. This question can be answered by ignoring for a moment the 
spatial sampling of the image that leads to pixelization. We imagine that the image intensity is 
a smooth function u{x) defined for x G [0, 1] x [0, 1]. In this case, V is clearly a two-dimensional 
smooth manifold. Even though the image is a function of a continuous argument a?, the patch u{x) 
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is formed by collecting m discrete samples on a lattice. Let A be the horizontal, or vertical, lattice 
resolution. For any pair of translation indices fc, I = — m/2, . . . , m/2, we consider the function that 
maps X to the ((/ + m/2) + m{k + m/2) + 1)^^ coordinate of the patch u{x)^ 



This map is a smooth function from [0, 1] x [0, 1] to R, since u is smooth. Therefore, each of the 
coordinates of u{x) is a smooth function of x, and the map that assigns a pixel to its patch 



is a smooth map from [0, 1] X [0, 1] to M"* . The non-discretized version of V is thus a two- 
dimensional manifold in M"^^. We note that Grimes and Donoho [17] have argued that a set of 
binary images obtained by moving a black object on a white background forms a sub-manifold that 
is not differentiable. The lack of differentiability hinges on the fact that binary images can only 
have step edges that are not differentiable. In fact, smoothing the images allows one to recover the 
differentiability of the sub-manifold. In this work, we assume that the image u has been blurred by 
the acquisition device, and therefore u{x) is a smooth function of x. 

We provide an illustration of the patch submanifold V in Fig. 1, where 5x5 patches, extracted 
from the 128 x 128 butterfly image (Fig. 1-left), form a smooth two-dimensional structure (Fig. 1- 
right). The cone-shaped surface encodes the content of the patches: as we move along the curve 
formed by the opening of the cone in Fig. 1-right, we explore patches (numbered 4, 5, 6 and 7) 
that contain high-contrast edges at different orientations. Each ruling (generator) of the cone is 
composed of all the patches that contain an edge with the same orientation. In effect, the orientation 
of the edge inside these patches is encoded by the position of the ruling on the opening curve. In 
addition, the white part of the patch is either always at the top, or always at the bottom (see e.g. 
4 — )► 3 1 and 6-^2-^1). Finally, as we slide from the opening of the cone to its tip, the contrast 
across the edge decreases (see e.g. 4 ^ 3 — > 1 and 6 2 -> 1). 

2.1. Denoising patch datasets 

We now consider the problem central to this paper: denoising a patch dataset that lies close to 
a manifold M. We review the two broad classes of approaches that have been proposed to solve 
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(4) 




(5) 




Figure 1: Manifold of 5 x 5 patches (right) extracted from the 128 x 128 butterfly image (left). Each dot represents a 
single patch u{x). A few patches are labeled in the image and on the patch manifold to assist with the interpretation 
of the manifold in terms of the geometric content of the patches. 

the more general problem of denoising a manifold- valued dataset: 

• local methods that construct a local projector, which varies from point to point, and project 
the noisy measurements on a local estimate of the tangent plane T^(a,)7W (see Fig. 2); 

• global methods that denoise the entire dataset. Linear and nonlinear methods can be used to 
reconstruct an estimate of the clean manifold. 

Local methods: denoising along the tangent plane. The concept of image patches is equiv- 
alent to the concept of time- delay coordinates in the context of the analysis of a dynamical system 
from the time series generated by an observable [18]. Several methods have been proposed to re- 
move the noise from time-delay embedded data [e.g., 19, and references therein]. These methods 
rely on the estimation of a local coordinate system formed by the tangent plane to the manifold 
at u{x) (see Fig. 2). The tangent plane is computed from the singular value decomposition of a 
neighborhood of u{x)] the first singular vectors provide a basis of the local tangent plane T^i^^^Ad 
(see Fig. 2), while the remaining vectors capture the noise [e.g., 20, and references therein]. The 
noise, which is dimensional, is squashed by projecting the noisy data onto the tangent plane 
(see Fig. 2). Similar ideas have been proposed recently in the image processing literature. The 
BM3D algorithm [5] uses the local patches around a noisy reference patch u(x) to construct a local 
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Figure 2: Local denoising along the tangent plane of the patch-set. The clean patch u{x) is corrupted by isotropic 
Gaussian noise. The noisy measurement u{x) is inside a ball centered around u{x). An estimate of u(x), is 
computed by projecting u{x) on the tangent plane to the manifold Tu(x)M- 

orthonormal basis. A multiscale version of this concept was recently proposed to denoise curves [21]. 
Elad and co-workers proposed the k-SVD algorithm [22] that assembles multiple singular vectors, 
computed locally, in order to construct a global set of functions to represent an image. 
Global methods: diffusion on the manifold. Instead of trying to approximate locally the 
manifold of patches by the tangent plane T^(^^^A4 around each patch u(x)^ global methods perform 
a comprehensive denoising on the entire manifold. The global denoising can be performed by 
applying a smoothing operator, defined by a kernel i7, to the noisy dataset and compute an estimate 
u(x) of u(x) from the noisy measurement 



In practice, the Gaussian kernel is usually used, and the geodesic distance on the manifold is 
replaced by the Euclidean distance. The recent nonlocal means algorithm [3] can be recast as a 
discrete implementation of (6). Szlam et al. [6] have explicitly modeled the manifold of patches with 
a graph and implemented various diffusions on the graph. Singer et al. [23] provided a stochastic 
interpretation of nonlocal means as a diffusion process. Similar ideas were proposed in [24, 7, 5]. 
Gilboa and Osher [4] proposed a continuous version of this idea, and defined PDE-based flows that 
depend on local neighborhoods. Early on, the computer graphics community has been implementing 
discrete versions of (6) to smooth and denoise surface meshes [e.g., 25, and references therein]. 
Specifically, Taubin [25] proposed to implement (6) using the diffusion 




(6) 



du{x, t) 
dt 



— XCu{x, t) 



(7) 



where £ is a discrete approximation to the Laplacian operator defined on the mesh. 



Global methods: from the diffusion operator to its eigenvectors. Instead of applying the 
diffusion kernel H on the graph of patches, one can compute the eigenvectors of the graph 
Laplacian [11, 12], and use them to expand the corresponding diffusion kernel H. The set of 
eigenvectors provides a spectral decomposition, similar to a global Fourier basis. If we keep only 
the first few eigenvectors of the expansion, the result is very similar to applying the diffusion 
operator H for a long period of time (see section 4.2). In his seminal work, Taubin [25] proposed to 
use the eigenvectors of the graph Laplacian to denoise meshes. Szlam et al. [6] proposed to perform 
some nonlinear approximation using the eigenvectors of the graph Laplacian. Peyre [8] combined 
thresholding and regularization to remove noise from images. 

In our problem, the diffusion operator H needs to be estimated from the noisy data. Conse- 
quently, we expect that the eigenvectors of the noisy operator H may be very different from the 
eigenvectors of the true operator H. We intend to study the perturbation of the eigenvectors, and 
propose a method to denoise jointly the eigenvectors and the image. 

3. The eigenvectors: a sparse code for images 

3.1. The graph of patches: the patch- graph 

In this paper, we think about a patch, u{xn)^ in several different ways. Originally, u{xn) is 
simply a block of an image. We also think about u{xn) as a point in . Finally, we also regard 
u{x) as a vertex of a graph. Throughout this work, we will use these three perspectives. In order 
to study the discrete structure formed by the patch-set (3), we connect patches to their nearest 
neighbors, and construct a graph that we call the patch-graph (see Fig. 3). 

Definition 3. The patch-graph^ F^ is a weighted graph defined as follows. 

1. The vertices of F are the N'^ patches u{xn)^ n = 1, . . . , N'^. 

2. Each vertex u{xn) is connected to its v nearest neighbors using the metric 



d{n,m) = \\u{xn) - u{xm)\\ + P\\xn - X 



•m 



(8) 
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The weight Wn,m along the edge {u{xn) , u{xm)} is given by 




(9) 
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Figure 3: The 3x3 patches u{x) are identified as vectors in R^. The patch-graph encodes the similarities between 
different regions of the images at a given scale. 

The parameter /3 > controls the influence of the penalty term that measures the distance between 
the patches in the image domain. The distance d{n, m) is small if the image intensity is similar at 
Xfi and u(^Xfi) ^ u^Xfy^)^ and and x^y^ are not too far from one another || Xfi XfYi '"^ 0. The 
parameter S controls the scaling of the similarity m) between u{xn) and u(xm) when defining 
the edge weight Wn^m- The particular form of the weight (9) ensures that Wn,m drops rapidly to 
zero as m) becomes larger than S. A very large S (i.e. Wn^m ^ 1 foi* all n, m) emphasizes the 
topology of the graph and promotes a very fast diffusion of the random walk through the patch-set. 
The other alternative: 6^0 accentuates the difference between the patches, but is very sensitive 
to noise. We now define the weight matrix, which fully characterizes the patch-graph. 

Definition 4. The weight matrix W is the N'^ x N'^ syuiTnetric uidtrix with entries — Wji . 

The degree matrix is the N'^ x N'^ diagonal matrix D with entries Dn^n — Ylm=i ^n,m- 

Finally, we define the normalized Laplacian matrix. 

Definition 5. The normalized Laplacian matrix L is the N'^ x N'^ symmetric matrix defined by 

L = I-D-^WD-^. (10) 
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We note that the sign of the Laplacian is the opposite of the sign of the Laplacian on a manifold 
M of dimension defined by 

where B{x^r) is the ball of radius r centered around and Yo\[B{x^r)) is the volume of the 
ball. Despite the fact that the discrete Laplacian has the wrong sign, it shares the structure of (11). 
Indeed, the entry n of the vector D~^WD~^u is an average of the function u defined on the graph, 
computed in a small neighborhood around the node u{xn)- Similarly, ^^^^^(^^ ^) u{y)dy 

computes the average of u within a ball centered around x G In both cases, the Laplacian 
measures the difference between u{x) and its local average. 

The matrix L is symmetric and positive definite, and has N'^ eigenvectors - , 0a/'2 with 

corresponding eigenvalues Ai = < Ai < • • • < A^ys . Each eigenvector (f)^ is a vector with A/"^ 
components, one for each vertex of the graph. Hence, we write 

iT 



0/c 



to emphasize the fact that we consider (pj^ to be a function sampled on the vertices of F. The 
eigenvectors 0i, . . . , 07v2 form an orthonormal basis for functions defined on the graph, where the 
inner product on the graph is defined by 

3.2. The eigenvectors encode the geometry of the image 

iT 



We can represent (f)^ = 



as an image: each pixel Xi is color-coded 



according to the value of (f)k{xi). Figure 4 shows the three vectors 02, 03 and 04 of the butterfly 
image (see Fig. 1-left). The first non trivial eigenvector 02 encodes the gradient of the image 
intensity: 02 takes large negative-^ values at pixels with large gradient (irrespective of its direction) 
and takes large positive values at pixels with small gradient. 03 and 04 clearly encode horizontal 
and vertical partial derivatives of u. The frequency content of each eigenvector 0/^; is loosely related 
to the eigenvalue A^. Unlike Fourier analysis, the basis functions 0^ have an 



^We note that the sign of (j)k is arbitrary. 
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Figure 4: From left to right: eigenvectors </)2, </)3, and cf)4. 



intrinsic scale given by the patch size. In this work we propose to use the eigenvectors <pk^k = 1, . . . 
as a basis to denoise the image intensity function u. 

We conclude this section with experimental evidence that indicates that the image u has a 
sparse representation in the basis {(pk}- We will use this property in the next section to remove 
the noise from the image. Figure 5 displays the relative residual energy, r^, after reconstructing 
the clown image (Fig 3-left) using the first K eigenvectors with the largest coordinates {u^ 4>k)-> 



\U 



TK 



\U\\ 



For comparison purposes, the same quantity is computed using a wavelet decomposition. We see 
that the basis functions 4>k give rise to a very fast decay of the residual energy, even faster than with 
the 9/7 wavelet basis. We performed similar comparisons with several other images (not shown) 
that lead us to conclude that we could use this basis to denoise natural images. 
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Figure 5: Relative residual energy after reconstructing the clown image using the i^-best eigenvectors (red), and 
using the K-hesi wavelets (green). 
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4. Denoising 

4.1. Denoising the patch-set using the noise-free eigenvectors 0^ 

We describe in this section a denoising algorithm that constructs an estimate of the patch-set 
from the knowledge of the noisy image u. The denoising of the patches is not performed locally on 
the patch-set, but rather relies on the global eigenvectors (pk- This allows us to denoise the entire 
patch-set as a whole. We assume for the moment that we have access to the eigenvectors <pk of 
the noise- free image. While this is clearly unrealistic, this allows us to explain the principle of the 
algorithm. In section 4.3 we study the effect of noise on the eigenvectors. Finally, in section 4.4, 
we describe a procedure to iteratively reconstruct the patch-set and the eigenvectors. 

The original image u has been corrupted by additive white Gaussian noise with variance cr^, 
and we measure at every pixel x the noisy image u{x)^ given by 

u{x) = u{x) + n{x). 

Let u{xn) be the patch centered at Xn associated with the noisy image 5, and let V be the patch- 
set formed by the collection of noisy patches. We propose to construct an estimate V of the clean 
patch-set V from the noisy patch-set V. Lastly, we will combine several denoised patches u{xm) 
to reconstruct an estimate u{xn) of the clean image at the pixel x^- 

Let ui{xn), ' ' ' ,Ur^2{xn) dcnotc the m? coordinates of the patch u{xn)- We define similarly, 
ui{xn)^ . . . ,5^2(x^) to be the m? coordinates of the patch G(x^). We make the following trivial 
observation: the coordinate uj is a real-valued function defined on the set of vertices 

y(r) — > M 
Uj '. Xfi I y Uj{xfi) 

In order to denoise simultaneously all the patches, we denoise each coordinate function Uj inde- 
pendently. This is achieved by expanding the function Uj into the basis formed by the 0^, and 
performing a nonlinear thresholding. The result is the denoised function Sj, given by 
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Figure 6: The image value u(i^j) (in red) is reconstructed using the coordinate of each of the rn? overlapping patches 
u{xi) that corresponds to the location Xn = {hj)- 

where the function hi performs a nonhnear thresholding of the coefficients (pk)- Iii practice, we 
use a hard thresholding. After denoising all the coordinate functions Si, . . . ,S^2, we have access 
to an estimate {u{xn)^n = 1, . . . , N^^ of the clean patch-set. 

Because the patch-set has a maximum overlap (i.e. we have a patch for each pixel), there are 
m? patches that overlap (in the image domain) with the pixel Xn (see Fig. 6). Let N{xn) be the 
ball for the /q norm of radius m/2 pixels centered at Xn-> 

N{Xn) = {xu \\xn - Xi% < m/2} . (13) 

The patches u{xi) that overlap with Xn all have their center xi in N{xn) (see Fig. 6). For 
each of the pixels xi G N{xn)-> there exists an index i^ such that the coordinate of the patch 
u{xi) corresponds to the pixel location Xn- In other words, [u{xi)\i^ is an estimate of the image 
intensity at the location Xn obtained from patch u{xi). We combine all these estimates and define 
the denoised image at pixel 

as the following; weighted sum 

^(^n) = ^ a{Xn,Xi)[u{xi)\i^, (14) 

Xi^N{Xn) 



a[Xn,xi)^— 11^ — l|2N • (15) 



where the exponential weight a(x^,x/) is given by 

exp {-\Xn - Xjf') 

This choice of weights favors the patch u(xn) centered at x^ and disregard exponentially fast 
the neighboring patches. We have also experimented with weights that discourage the mixing of 
patches with very different variances (results not shown). 

Connection to translation-invariant denoising. We note that we can interpret the jth coor- 
dinate function uj as the original image u shifted by (p, g), where j = m(p + m/2) + g + m/2 + 1. 
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For instance, if j = 1, then the two-dimensional shift is (— m/2, — m/2) and if j = m^, the shift 
is (m/2,/m2). The reconstruction formula (14) can then be interpreted as a translation-invariant 
denoising procedure [26] of the form, 

u{x) = Average < Unshift Denoise Shift u{x) > . 

Me[-m,m\x[-m,m] I M (m) J 

4-2. Geometric interpretation of the eigenvectors 

We explain in this section the connection between our approach and denoising methods based 
on defining a diffusion process on the graph [6, 7]. The solution to the denoising problem is related 
to a diffusion process defined on the patch-set. Indeed, let us consider the image v(x^t) solution 
to the diffusion defined on the patch-graph (notice that we use the wrong sign for equation (7) as 
explained in section 3.1) defined by, 

^ = -Mx,t), (16) 

where the initial condition is provided by the noisy image u{x)^ 

v{x,0) = u{x). (17) 
We consider the heat kernel formed by the N'^ x A^^ matrix 

Ht{n, m) = Y^ e-^^^ct)k{xn)ct)k{xm)^ (18) 
and we write the solution to the diffusion as 

Ar2 7V2 

v{xn, ^) = XI Ht{n, m)u{m) = ^ e"^^* (2, 0^) (t>k{xn). (19) 

171=1 k=l 

We notice the similarity between the expansions (12) and (19). As t becomes large, only a small 
number of terms in (19) will be non zero. The attenuation of (2, 0/c) is controlled by the size 
of the corresponding eigenvalue Xk. This is in contrast to (12), where the coefficients used to 
reconstruct the denoised signal are chosen based on the magnitude of (5, The similarity 

between the two expansions allows us to understand what is the information encoded by each (p^. 
In (19), each function e~^^^ (t)k{x) is a stationary solution to the diffusion (16). The amplitude of 
this stationary solution decreases exponentially fast with time, but the geometry, encoded by (j)]^ 
remains unchanged. The eigenvectors 0^ therefore encode the geometric features present in the 
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image: edges, texture, etc. as suggested by Fig. 4. Finally, we note that if the graph of patches is 
replaced with a regular lattice formed by the image sampling grid, and if each patch is reduced to a 
single pixel, then the diffusion (16) models one of the standard diffusion-based denoising algorithms 
[27]. 

4.3. Estimating the eigenvectors 4>k from the noisy images 

We described in the previous section a denoising procedure that relies on the knowledge of the 
eigenvectors 0^ associated with the clean patch-set. Of course, if we had access to the clean patch- 
set, then we would have the clean image. While our approach appears to be circular, we claim 
that it is possible to bootstrap an estimate of the (f)^ and iteratively improve this estimate. This 
approach relies on the stability of the low- frequency eigenspace of the graph Laplacian L. In the 
next section we study experimentally the perturbation of the eigenvectors 0^ when a significant 
amount of noise is added to an image. We demonstrate that changes in the topology of the graph, 
and not changes in the weights, is the source of the perturbation of the (pk- This result leads to an 
iterative method that recovers the clean eigenvectors (j)^ that is described in section 4.4. 

4.3.1. The eigenvectors of the perturbed Laplacian 

Each patch u{xn) of the clean image is corrupted by a multivariate Gaussian noise n{xn) with 
diagonal covariance matrix, 

u{Xn) = u{Xn) + n{Xn). (20) 

The distance between any two noisy patches u{xn) and u{xrm) is given by 

\\u{xn)-u{xm)\\^ = ||u(x^) - u(x^) f + 2(u(x^) - u(x^) , n(x^) - n(x^)) + \\n{xn) - n{xm)\\^ ' 



(21) 

with the corresponding expected value, 

^\u{Xn) - u{Xm)f = \\u{Xn) - u{Xm)f + 2aW. (22) 



The weights associated with the noisy image are given by 

While e-|l^(«^n)-n(a.^)||2/5' < the term e-2{^(a^n)-tx(a.m),n(a.n)-n(a.m)>/(52 greater than 1, 

thereby increasing Wn,m- Overall, the presence of noise has two effects: 
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1. the topology of the graph is modified since points that were previously neighbors may become 
far away and vice- versa (caused by changes in the distance (22)), 

2. the weights along the edges are modified according to (23). 

Let r be the graph associated with the noisy patch-set V. V and V have the same vertices, 
V{V) = V{V)^ and the same total number of edges (since both are based a nearest neighbor 
topology). The weight matrix W is perturbed according to (23). We define W and L to be the 
weight and Laplacian matrices computed from the noisy image, respectively. We expect that the 
eigenvectors (j)^ of L will be different from the (f)]^- As explained in the next section, the low- 
frequency eigenvectors remain stable, but the higher frequency eigenvectors rapidly degrade. As a 
result, 0fc cannot be used to reconstruct a clean image. 

Theoretical results that provide upper bounds on the angle between an eigenspace of L and 
the corresponding eigenspace of L exist [13]. Unfortunately, as noted in [14], these bounds usually 
overestimate the actual perturbation of the eigenvectors. Furthermore, these results depend on 
the separation of the eigenvalues, a quantity that is difficult to predict in our problem. In fact 
the separation can be very small leading to very large bounds. Another limitations of the usual 
theoretical bounds is that the angle between two invariant subspaces cannot be readily translated 
in terms of the ability of the perturbed eigenvectors (pk to approximate the original image using a 
small number of terms. We propose therefore to study experimentally the effect of the noise on the 
ability of 0^ to encode geometric information about the image. We employ standard methodology 
used in vision research. Specifically, let J^(f)k be the Fourier transform of 0/^, we study the energy 
distribution of J^(f)k using polar coordinates (p, 6) in the Fourier plane and we compute the average 
energy of J^4>k over all orientation 6 for a given radius p. 

Our experiment was performed as follows. We add white Gaussian noise (cr = 40) to the clean 
clown image (see Fig. 7), and we compute the eigenvectors (pk of the graph Laplacian. This image 
provides a good balance between texture and smooth regions; we performed similar experiences with 
other images, and obtained results that were quantitatively similar. The perturbed eigenvectors 
02, 032? 01285 ^ud 0256 ^1"^ shown in the bottom row of Fig. 8. The eigenvectors of the clean image 
are shown in the top row of Fig. 8 for comparison purposes. We can make the following visual 
observations, which will be confirmed by a quantitative analysis in the next paragraph. For k < 32, 
the eigenvectors <pk appear to be reasonably similar to the original eigenvectors (pk- For k > 128, 
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original noisy 




Figure 7: Original 128 x 128 "clown" image (left); clown image with white Gaussian noise, a = 40 (right). 



<p2 032 0128 0256 




02 032 0128 0256 

Figure 8: Top row: eigenvectors of the clean image. Bottom row: eigenvectors of the noisy image. 

4>k appear to mostly capture the random noise instead of the edges and texture present in the 
original image. For instance, the texture on the wallpaper and on the shirt of the clown, which is 
present in 0256, is replaced by noise in 0256 (see Fig. 8). 

A quantitative analysis of the perturbation of the 0/^. We provide here a quantitative 
comparison between 0/^ and (pk- Our evaluation is based on the comparison between the energy 
distribution of the Fourier transforms Tcj^k and J^0/c- Following a practice standard in vision 
research, we analyse Tc/^k as a function of the radial and angular frequencies, (p^O). Because we 
assume that the visual features (edges, textures, etc) can have any random orientation, we integrate 
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\J^(t>k\^ over all possible orientation 9 and compute the energy of J^cj^k within a thin ring i?/, 

^fc(0= / \F^k\\p.0)dpde, (24) 
JRi 

where the ring Ri is defined by 

Ri^{{p^e) G [/Ap, (/ + 1) A^] X [0,2^]}, Z = 0,...,L-1. 

The width of each ring (Ap) depends on the sampling frequency and the size of the image. In the 
experiments we use a total of L = 32 rings, and Ap = 2 for an image of size 128 x 128. In other 
words, there were two discrete frequencies within each ring i?/. We do not study the perturbation 
of each eigenvector individually. Instead, we group the eigenvectors (pk according to a dyadic scale 
of the index k: the scale i corresponds to the eigenvectors 

, 02^+1 } ^ i = 0, .... 

This dyadic repacking of the eigenvectors was motivated by the experimental observation that the 
eigenvectors 0^/^ with the same dyadic scale i had a similar Fourier transform. Finally, we compute 
at each scale i = 0, . . . the total energy contribution, integrated inside the ring i?/, from the group 
of eigenvectors at that scale i. 

Figure 9 shows the energy distribution £'^{l) as a function of the radial frequency index /, for the 
scales i = 0, . . . 7. For each plot, the x-axis is the index I of the ring Ri, and therefore provides a 
quantization of the radial frequency. The y-axis is the total energy contribution, measured inside 
the ring i?/, from the group of eigenvectors at scale i. We plot the energy £'^{l) of the perturbed 
eigenvectors (in blue), as well as the energy £'^{l) of the clean eigenvectors (in red). 

For the first three scales i = 0, 1, 2, corresponding to = 2, . . . , 8, there is hardly any difference 
between the energy distribution of the clean and perturbed eigenvectors (see Fig. 9). Starting at 
scale i = 3, the energy of the perturbed eigenvectors becomes different from the energy of the clean 
eigenvectors: has a flatter distribution with more energy in the higher frequencies. The energy 
leak into the high radial frequencies is created by the noise and confirms the visual impression 
triggered by 0256 Fig- 8. The high-index eigenvectors (f)^ are trying to capture the noise present 
in the image u. As the scale i further increases, the departure of from becomes even more 
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Figure 9: Energy S^{1) and S^{1) of the eigenvectors (pk and </)fc, respectively, as a function of the radial frequency 
index / for the "clown" image (i = 0, . . . 7 from top to bottom, and left to right). The cf^k are computed using the 
noisy image. 



dramatic. At scale i = 7, the radial energy distribution of the eigenvectors (/)i29, • • • , </>256 is almost 
flat, indicating that the (pk have adapted to the noise. These eigenvectors are obviously no longer 
suitable for denoising. 



Which perturbation matters most: the topology of the graph, or the edge weights? 

As demonstrated in some recent studies [28, 29] topological changes of a graph (created by the 
addition or removal of edges) can significantly perturb the spectrum and the eigenvectors of the 
graph Laplacian. Inspired by these studies, we analyze the effect of the modifications of the topology 
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of the graph on the eigenvectors 0^. 

Precisely, we verify experimentally the following result: if the non-zero entries in the matrix W 
(9) are kept at the same locations (we keep the graph topology), and if we randomly perturb the 
weights according to (23), then the perturbed eigenvectors have almost the same energy distribution 
in the Fourier domain as the eigenvectors computed from the original weight matrix. In the previous 
paragraph we observed that random fluctuations of the weights combined with perturbations of the 
graph topology (caused by the random changes of the local distance, and the corresponding neigh- 
bors) significantly perturbed the eigenvectors of L. Combining these two experiments, we conclude 
that the topology of the graph is the most important factor, since, if preserved, it guarantees the 
stability of the eigenvectors. 

Let us now describe the experiment. We first build a patch-graph and a weight matrix W 
based on the clean image. We now add a random realization of white Gaussian noise {a = 40), 
and compute the new weights Wn,m according to (9), where the neighbors are now defined by the 
clean patch-graph. We then update all the non-zero entries of the matrix W using the perturbed 
weights Wn,m- This constitutes the perturbed matrix W. We define L accordingly. 

We note that W is not equal to the matrix described in the previous paragraph. Indeed, 
the nonzero entries of W correspond to the edges of the graph defined by the perturbed distance 
dn^m- In other words, Wn^m — ^n,m oiilj if u{xn) and u{xm) are neighbors, and u{xn) and u{xm) 
are also neighbors. 

Let 0^ be the eigenvectors of the perturbed matrix L. Figure 10 displays 4>2-, ^325^1285 
0256- Despite the fact that the entries of the matrix L have been perturbed by adding a large 
amount of noise to the image, the eigenvectors of the resulting matrix, L, appear visually very 
similar to the original eigenvectors of L. This observation is confirmed in Fig. 11, which shows 
the energy distribution £'^{l) as a function of the radial frequency for the two sets of eigenvectors: 
(f)k and Remarkably, the distribution of energy of (j)^ and 0^, measured across the different 
radial frequencies in the Fourier plane, are very similar. The eigenvectors (p^ do not suffer from 
any significant perturbation despite very significant changes in the image patches. 

This experimental result is related to similar properties in the context of spectral geometry. If 
we replace the graph F by a manifold, and we consider that the weight matrix encodes the metric, 
then our experimental results are related to the more general phenomenon that involves the stability 
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02 032 0128 0256 

Figure 10: The topology of the graph that is used to compute the is determined from the clean image. 
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Figure 11: Energy 8^ {I) and of the eigenvectors (red) and (pk (blue) as a function of the radial frequency 

index / for the "clown" image (i = 0, . . . 7 from top to bottom, and left to right). The topology of the graph that is 
used to compute the 4>^ is determined from the clean image. 
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of the eigenspaces of the Laplace-Beltrami operator when the metric on the manifold is perturbed 
continuously [15, 16]. In our case, the continuous changes of the metric correspond to changes in 
the weights Wn^m- 

In practice, we do not have access to and therefore we cannot construct W. However, the 
results of this section suggest that we should try and reconstruct an estimate of the graph topology, 
associated with W. We combine the following two results to construct an estimate of W. 

1. For most images, many of the entries of the weight matrix W can be estimated from a lowpass 
version of the image. In other words, 

\\u{Xn) - u{Xm)\\ - \\hu{Xn) - hu{Xm)\l 

where /i is a lowpass filter, and hu is a patch from the image hu. Indeed, unless we are 
studying patches that contain only high frequencies, most of the energy within the patch 
comes from the low frequencies. 

2. The eigenvectors from the very first scales, i — 0, . . . , 4, are very stable even when a large 
amount of noise is added to the image. We can therefore use these eigenvectors to estimate a 
lowpass version of the original image. 

The strategy for denoising is now clear. We compute the weight matrix W and calculate the 
eigenvectors associated with the first scales i = 0, 1,2,3,4. We use these eigenvectors to compute 
a coarse estimate of the lowpass part of the image, u^^\ using the procedure described in section 
4.1. Because too few eigenvectors are used in this reconstruction, we add back a scaled version of 
the residual error (the difference between the noisy image u and the denoised image vl^^^) to obtain 
the intermediate image S^^^ defined by 

= + ^{u - = (1 - 7)?x(^) + 7S, (25) 

where 7 G (0, 1). at this point the image iP'^ has been sufficiently denoised, and we construct a 
new graph. We note that we have experimented with a different procedure to compute the image 
u^^\ using standard lowpass filtering, or wavelet denoising. Using these alternate approaches, we 
obtain similar results, albeit of lower quality (results not shown). 

We finally compute the eigenvectors 0/^ of L^^^ associated to the image u^'^\ Figure 12 displays 
the eigenvectors 02, 032, 0128 5 ^iid 0256- These eigenvectors are visually similar to those 

21 





02 032 0128 0256 

Figure 12: The eigenvectors (pk are computed using a second graph that is built from the denoised image u^'^\ 
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Figure 13: Energy S^{1) and S^{1) of the reconstructed eigenvectors (pk (red) and <pk (blue) as a function of the radial 
frequency index / for the "clown" image (i = 0, . . . 7 from top to bottom, and left to right). The eigenvectors <pk are 
computed using a second graph that is built from the denoised image u^'^\ 
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computed using the graph of the clean image (see Fig. 10). Figure 13 displays the distribution of 
the energy of the eigenvectors 4>k as a function of the radial frequency. Even at large scale i = 6 
and 7, the distribution of the energy of the cj)k across the radial frequency almost coincides with 
the distribution of the energy of the (pk- 

Szlam et al. [6] proposed a related approach to bootstrap a semi-supervised classifier: a new 
graph is constructed after one iteration of the classifier to improve the classification performance. 

4.4- Iterative denoising 

We propose an iterative procedure to denoise an image based on the results of the previous 
section. In practice we use two passes of denoising. The algorithm is described in Fig. 14. In the 
next section we study in details the two most important parameters of the algorithm: the patch 
sizes (mi and m2) and the numbers of eigenvectors {Ki and K2) that are used during the two 
stages of the algorithm. 

4.5. Optimization of the parameters of the algorithm 
4.5.1. The patch size 

While the eigenvectors 0^ provide a global basis on the patch-set, the patch size introduces 
a notion of local scale in the image. We notice visually that as the patch size increases, the 
eigenvectors become less crisp and more blurred. As explained in [23], [9] when the patch size is 
larger, patches are at a larger distance of one another in the patch-set. Consequently, a larger 
patch size prevents patches that are initially different from accidentally becoming neighbors after 
the image is corrupted by noise, since their mutual distance is relatively less perturbed by the 
noise. Because the eigenvectors <pk are sensitive to the geometry of the patch-set, larger patches 
help minimize the perturbation of the eigenvectors. We note that, as patches become larger, the 
eigenvectors become more blurred, therefore requiring more eigenvectors to describe the texture. 
Finally, we need to be aware of the following practical limitation: if the patch size is large, and the 
image size is relatively small, then we have few patches to estimate the geometry of the patch-set; 
a problem aggravated by the fact that the patch-set lives in a large ambient space in that case. 

We have compared the effect of the patch size on the quality of the denoised image S^^^ after 
the first pass (stage 1) of the denoising algorithm (see Fig. 14). We define the signal to noise ratio 
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Algorithm: Two-pass denoising 



Input: noisy image u; 

• number of nearest neighbors patch sizes mi,m2; scale parameters 61^62] 

• number of eigenvectors Ki, K2] 7 G [0, 1] 

Stage 1 

1. build the graph from the image u and compute W and D; 

2. compute the first Ki eigenvectors, <pk, = 1, . . . of D~2WD~2 

3. construct the nonlinear estimator u^^\x) of the patch-set using (12) 

4. reconstruct the image u^^"^ by averaging the patches u^^\x) using (14) 

5. compute the denoised image S^^^ = (1 — 7)^^-'^) + 

Stage 2 

1. build a second graph from the image 2^^^; compute Vl^(^) and D^^); 

1 

2. compute the first K2 eigenvectors, cl)k, A: = 1, . . . , K2 of 

3. construct the nonlinear estimator u^^\x) of the patch-set using (12) 

4. reconstruct the image S^^^ by averaging the patches u^^\x) using (14) 

Output: the denoised image u^^\ 



Figure 14: The two-stage denoising algorithm. 

(SNR) by 

SNR= „ " (26) 
\\u — u\\ 

Figure 15 displays the signal-to-noise ratio as a function of the number of eigenvectors K used to 
reconstruct u^'^\ for different patch sizes m, and at different noise levels {a = 20,40, and 60). For 
low noise levels (a = 20) the SNR is larger for smaller patches. For moderate noise levels, the 
maximum SNR achieved with different patch size is approximately the same (see Fig. 15 bottom- 
right). Nevertheless, the visual quality of the denoised images - at the same SNR - are quite 
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Figure 15: Top to bottom, left to right: SNR as a function of the number of eigenvectors K after a one pass 
reconstruction, for increasing noise levels. Each curve corresponds to a different patch size. Bottom right: optimal 
reconstruction error as a function of the patch size m at different noise levels. 




Figure 16: Images with the same SNR denoised with patch size m = 3, 5, 7, 9 (from left to right). 

different. Figure 16 shows four denoised images (a = 40) with very similar SNR, for m = 3, 5, 7 and 
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9. With m = 3, the artifact caused by the noise can still be seen in the denoised image. However, 
with larger m, the denoised image is much smoother, although some of the details in the original 
image are lost. The choice of the patch size m is in general based on the noise level a. In practice, 
we choose m = 3 for a = 20, m = 5 for a = 40, and m = 7 for a = 60. 

4-5.2. The number of eigenvectors K 

At both stages of the algorithm, the first K eigenvectors are used to compute a denoised image. 
In general, the eigenvectors degrade faster as the variance of the noise level increases. We explain in 
the following why, using a very small number of eigenvectors, we can still reconstruct the structures 
of the patch-set that contain patches extracted from smooth regions of the image. 
The local dimensionality of the patch-set. In this work we consider images that contain a 
mixture of smooth regions and texture. Therefore, we only need a very small number of eigenvectors 
to capture the patches extracted from smooth regions of the image. Indeed, as explained by Taylor 
and Meyer [9], patches extracted from smooth regions of the image align themselves along smooth 
low-dimensional structures. Indeed, for smooth patches, spatial proximity (in the image) implies 

2 

proximity in . Conversely, patches that contain large gradients and texture are usually scattered 
across M"^^ [9]. Other researchers have made similar experimental observations. For instance. 
Chandler and Field [30], Zontak and Irani [10] observed that patches can be classified into two 
classes according to their mutual distances. One class is composed of patches extracted from smooth 
regions, and the other class encompasses patches extracted from regions that contain texture, fast 
oscillation, large gradient, etc. These authors observe that "smooth patches" are at a small distance 
(in the patch-set) of one another, and are also at a close spatial distance (in pixels). They also note 
that proximity between "texture patches" is not dictated by spatial proximity (in pixels). Texture 

2 

patches are scattered in because they are at a large distance of one another: Zontak and Irani 
[10] estimate that the distance to the nearest neighbor of a patch grows exponentially with the 
gradient within that patch. Finally, Taylor and Meyer [9] have shown that the local dimensionality 
of the patch-set is much lower around smooth patches, and much higher around texture patches. 

Consequently, the smooth low-dimensional structures formed by the smooth patches can be 
accurately represented with a very small number of eigenvectors. Because low-index eigenvectors 
are very stable, we are able to reconstruct the smooth patches using the first few eigenvectors (pk of 
L. This partial reconstruction can then be used to estimate a new patch graph. The reconstruction 
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Figure 17: Reconstruction error of the clown image, after a single pass of denoising, as a function of the number of 
eigenvectors at various noise levels. The patch size m was 5 except for a = 60 where it was 7. 

of the texture patches requires many more eigenvectors because that part of the patch-set has a 
much higher local dimensionality. Using the second graph, one can then compute more stable 
eigenvectors, and reconstruct the texture patches. 

Figure 17 displays the SNR of u^^^ after a one-stage reconstruction of the clown image as a 
function of the number K of eigenvectors, for several noise levels. The patch size was m = 5, 
except for a = 60 where it was m — 1 . As we increase the number of eigenvectors the SNR 
initially increases. Except at very low noise level (a < 20), the SNR decreases after it reaches 
its maximum. Indeed, at moderate and high noise levels, the high index eigenmodes adapt to the 
incoherent structure of the noise and become useless for the reconstruction. This is why we need 
a second pass in the reconstruction. Figure 17 can help us determine the numbers Ki and K2 of 
eigenvectors needed for the first and second denoising passes. 

5. Experiments 

We implemented the two-stage denoising algorithm (see Fig. 14) and evaluated its performance 
against several denoising gold-standards. For all experiments, we kept Ki — 35 and K2 — 275 
eigenvectors during the first and second passes of the algorithm, respectively. The patch sizes for 
the first pass were mi = 7 for a = 40, and mi = 9 for a = 60, respectively. We used m2 = 5 for all 
noise levels in the second pass of the algorithm. 

The evaluation was performed using five images that contain a mixture of smooth regions, 
periodic and aperiodic texture (see Fig. 18). White Gaussian noise was added to the images. We 
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Figure 18: Top: clean images "clown", "lena", "roof", "barbara", and "mandrill". White Gaussian noise with cr = 40 
(middle) and a = 60 (bottom), is added to the images. 



used two noise levels: moderate a = 40 (Fig. 18 middle row), and high a = 60 (Fig. 18 bottom 
row). We compared our approach to the following three denoising algorithms: 

1. translation invariant soft wavelet thresholding (TISWT) method [26]. We used the imple- 
mentation provided by Wavelab (http://www-stat.stanford.edu/~wavelab/). 

2. nonlocal means [3]. We used the implementation provided by Jose Vicente Manjon Herrera 
and Antoni Buades (http://personales.upv.es/jmanjon/). 

3. k-SVD algorithm [22]. We used the implementation provided by Ron Rubinstein (http: 
//www. cs . technion. ac . il/~ronrubin). 

Table 1 displays the mean squared error between the reconstructed image and the original images, 

1 ^ 

^5^|S(^)(z,j)-^(^,j)P. 

^,J=l 

At all noise levels, and for all images our algorithm outperformed the other algorithms in terms of 
mean squared error. In terms of visual quality, the wavelet denoising yielded consistently the worst 
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noise level 


Laplacian 


k-SVD 


NL-means 
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40 


216 
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359 
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354 


509 
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236 
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329 
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298 


413 


382 


529 




40 


196 


318 


221 


486 
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310 


542 


395 


767 




40 


116 


160 


168 


241 


barbara 


60 


191 


300 


290 


436 




40 


314 


366 


357 


418 


mandrill 


60 


442 


542 


480 


621 



Table 1: Mean squared error (smaller is better) for the four denoising algorithms. 



reconstruction. Missing wavelet coefficients are very noticeable, even at moderate noise level. K- 
SVD and nonlocal means yielded similar results. The nonlocal means estimate was always more 
noisy than the k-SVD estimate, which was often too smooth. In comparison, our approach could 
restore the smooth regions and the texture (see e.g. the mandrill image at a = 40 and a = 60), 
even at high noise level. 

6. Discussion 

We proposed a two-stage algorithm to estimate a denoised set of patches from a noisy image. 
The algorithm relies on the following two observations: (1) the low-index eigenvectors of the diffu- 
sion, or Laplacian, operators are very robust to random perturbations of the weights and random 
changes in the connections of the graph; and (2) patches extracted from smooth regions of the 
image are organized along smooth low-dimensional structures of the patch-set, and therefore can 
be reconstructed with few eigenvectors. Experiments demonstrate that our denoising algorithm 
outperforms the denoising gold-standards. This work raises several questions that we address in 
the following. 
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6.1. Fast computation of the eigenvectors 

A key tenet of this work is our ability to compute the eigenvectors of the sparse matrix L, which 
has size N'^ x N'^ but with only u non zero entries on each row. For the experiments we used the 
restarted Arnoldi method for sparse matrices implemented by the Mat lab function eigs to solve 
the eigenvalue problem. 

There are several options for further speeding up the computation of the eigenvectors. Saito 
[31] proposes to modify the eigenvalue problem, and compute the eigenvectors via the integral 
operator that commutes with the Laplacian. This approach leads to fast algorithms. The recent 
work of Kushnir et al. [32] indicates that multigrid methods yield an immediate improvement over 
Krylov subspace projection methods (e.g. Arnoldi's method). Another appealing approach involves 
fast randomized methods [33]. We note that the application of these methods still necessitates the 
computation of the entire matrix W. We believe that more efficient methods that randomly sample 
the patch-set should be developed. While the sampling density should clearly guarantee that short 
wavelength (high-index) eigenvectors are adequately sampled, we can probably estimate the low- 
index eigenvectors with a very coarse sampling. Similar ideas have been recently proposed for the 
reconstruction of manifolds [34]. 

6.2. Extension 

This work opens the door to a new generation of image processing algorithms that use the 
eigenvectors of the graph Laplacian to filter the patch-set. We note that some of the most successful 
inpainting and super-resolution algorithms already operate locally on the patch-set [e.g., 35, and 
references therein]. These algorithms currently do not take advantage of the existence of global basis 
functions to represent the patch-set. Additional extensions include the construction of patch-sets 
from large collections of images. Lee and Mumford [36] demonstrated that high-contrast patches 
extracted from optical images were organized around 2-dimensional smooth (C^) sub-manifold. 
Simple models of synthetic images were constructed in [8] and were shown to lie close to low- 
dimensional manifolds. We note that the set of all m x m image patches includes the sub-manifold 
of mxm patches constructed from a single image, which is the sub-manifold studied in this paper. 

6.3. Open questions 

This work provides experimental evidence of the stability of the low-indices eigenvectors of the 
Laplacian L defined on a graph. Similar results exist in the literature in the context of the stability 
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Laplacian k-SVD NL-means TISWT 




Laplacian k-SVD NL-means TISWT 

Figure 19: Noise level a = 40. From left to right: our two-stage approach, k-SVD, nonlocal means, translation 
invariant wavelet denoising. 
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Laplacian k-SVD NL-means TISWT 




Laplacian k-SVD NL-means TISWT 

Figure 20: Noise level a = 60. From left to right: our two-stage approach, k-SVD, nonlocal means, translation- 
invariant wavelet denoising. 
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of the eigenvectors of the Laplace-Beltrami operator when the domain is perturbed [37], or when 
the metric is perturbed [15, 16]. To the best of our knowledge, there appears to be httle work on 
more precise analysis of the perturbation of the eigenspaces of the graph Laplacian. Yan et al. [14] 
acknowledge that standard bounds (e.g. [13]) overestimate the perturbations of the first eigenvector 
of the matrix L. 
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